home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / kaos / kaos_pat.0 / fp_get_evinfo.c next >
C/C++ Source or Header  |  1991-09-30  |  2KB  |  70 lines
  1. /* Order eigenvalues in decreasing order in magnitude */
  2. #include "../modellib/class_kaos_def.h"
  3.  
  4. #define DELTAX 1.e-7
  5. fp_get_evinfo(type,eval,evec,x,r,n)
  6. int *type,r,n;
  7. double **eval,**evec,*x;
  8. {
  9.     int i,j,*iv1,*ivector(),ierr;
  10.     double fabs(),sqrt(),**alpha,*beta,*x2,*dvector(),**dmatrix();
  11.     double *eval_r,*eval_i,**evectors,*fv1,**a;
  12.     extern int model,mapping_on,fderiv_on;
  13.  
  14.     x2 = dvector(0,n-1);
  15.     beta = dvector(0,n-1);
  16.     alpha = dmatrix(0,n-1,0,n-1);
  17.     eval_r = dvector(1,n);
  18.     eval_i = dvector(1,n);
  19.     evectors = dmatrix(1,n,1,n);
  20.     a = dmatrix(1,n,1,n);
  21.     iv1=ivector(1,10000);
  22.     fv1 = dvector(1,10000);
  23.  
  24.     for(i=0;i<n;i++) x2[i] = x[i]+DELTAX;
  25.     if(mapping_on){
  26.         if(fderiv_on){
  27.             usrfun(x,alpha,beta,r,n);
  28.         }
  29.         else{
  30.             usrfun2(x,x2,alpha,beta,r,n);
  31.         }
  32.         /* need this since alpha = D(f^r(x)-x) = Df^r(x) - I */
  33.         for(i=0;i<n;i++) alpha[i][i] += 1.;
  34.     }
  35.     else {
  36.         usrfun3(x,x2,alpha,beta,r,n);
  37.     }
  38.     if(n==1){
  39.         eval[0][0] = alpha[0][0];
  40.         eval[0][1] = 0;
  41.         evec[0][0] = 1;
  42.     }
  43.     else {
  44.         for (i=0;i<n;i++)
  45.           for (j=0;j<n;j++)
  46.              a[i+1][j+1] = alpha[i][j];
  47.         ierr=rg(n,n,a,eval_r,eval_i,1,evectors,iv1,fv1);
  48.         for(i=0;i<n;i++){
  49.             eval[i][0] = eval_r[i+1];
  50.             eval[i][1] = eval_i[i+1];
  51.             for(j=0;j<n;j++)
  52.              /*   evec[i][j] = evectors[i+1][j+1];   mrm  093091 */
  53.                 evec[j][i] = evectors[i+1][j+1];  
  54.         }
  55.     }
  56.  
  57.     /* get a type of a periodic orbit from eigenvalues */
  58.     fp_get_type(type,eval,n);
  59.  
  60.     (void) free_dvector(eval_r,1,n);
  61.     (void) free_dvector(eval_i,1,n);
  62.     (void) free_dvector(fv1,1,n);
  63.     (void) free_imatrix(iv1,1,n);
  64.     (void) free_dmatrix(evectors,1,n,1,n); 
  65.     (void) free_dmatrix(a,1,n,1,n); 
  66.     (void) free_dvector(beta,0,n-1);
  67.     (void) free_dvector(x2,0,n-1);
  68.     (void) free_dmatrix(alpha,0,n-1,0,n-1);
  69. }
  70.